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ABSTRACT 


dtic  quality  INSPECTED  3 


The  problem  of  linearly  constrained  least  squares  has  many  applications  in  signal  processing.  In  this  paper, 
we  present  a  perturbation  analysis  of  a  linearly  constrained  least  squares  algorithm  for  adaptive  beamforming.  The 
perturbation  bounds  for  the  solution  as  well  as  for  the  mfest  residual  element  are  derived.  We  also  propose  an  error 
estimation  scheme  for  the  residual  element,  which  can  be  incorporated  into  a  systolic  array  implementation  of  the 
algorithm. 

1.  INTRODUCTION 

The  least  squares  problem  with  linear  equality  constraints  has  important  applications  in  signal  processing, 
e.g.,  adaptive  beamforming.  To  solve  this  problem.  McWhirter  and  Shepherd  [5]  proposed  a  systolic  algorithm  and 
architecture.  In  this  paper,  we  present  a  perturbation  analysis  of  the  problem  and  propose  an  error  estimation 
scheme  for  the  McWhirter-Shepherd  (MS)  algorithm  [5].  This  paper  is  organized  as  follows.  The  least  squares 
problem  is  defined  in  Section  2  and  error  bounds  are  derived  in  Section  3.  An  error  estimation  algorithm  is  given 
in  Section  4,  and  in  Section  5  a  numerical  example  is  presented  to  illustrate  how  well  our  new  algorithm  works. 

2.  PROBLEM  DEFINITION 

Given  an  n  x  q  complex  data  matrix  A(n),  the  least  squares  problem  with  linear  equality  constraints  is  to  find 
a  (/-element  complex  vector  u'(n)  such  that 


|A'(n)u.'(n)||  =  min 


(2.1a) 


subject  to  the  linear  constraints 


Sw(n)  =  b. 


where  S  is  a  k  x  q  ( k  <  q)  complex  matrix  and  6  is  a  I’-element  complex  vector.  Throughout  this  paper,  we  use 
the  2-norin; 

II  II  =  11  lb- 

In  signal  processing,  new  data  arrives  continuously.  Define  the  data  matrix  A"(n)  recursively  by 


i.e.,  the  nth  row  i(n)T  represents  a  snapshot  at  time  n.  Our  goal  is  to  compute  the  n-th  residual  element 

rn  =  x(»)Tn(n). 
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Is  the  solution  vector  tr(n)  unique?  Define  &  (k  +  n)  x  q  matrix  S\(n)  by 


S x(n)  = 


We  assume  that  h  +  n  >  q.  The  solution  is  unique  if  and  only  if  the  matrix  S\ { n )  has  full  column  rank:  that  is. 
the  overdetermined  matrix  equation 

Sx  (n  )«’(«)  —  0  (2.3) 


has  a  unique  sohition  w(n)  =  0. 

Next,  we  wish  to  transform  (2.1)  into  a  familiar  unconstrained  problem:  see  [3]  and  [4j.  Let 


p  =  q  ~  k 


and  partition  the  matrix  .9  as 


9  =  (  9,  So  ) . 


where  9'i  is  lc  x  lc  and  So  is  h  x  p.  For  simplicity,  we  assume  that  9i  is  nonsingular  and  upper  triangular:  for 
example,  9'i  may  be  the  result  of  an  initial  QR  decomposition  of  9.  Accordingly,  we  also  partition  A  (n)  as 

,V(n)  =  (  A’i(n)  A’,(n) ), 


so  that  A'i  is  n  x  lc  and  A’o  is  n  x  p.  Then  (2.3)  becomes 


9i 

A’i(n) 


9,  \ 

A  2<  n )  ) 


w{n)  =  0, 


which  is  equivalent  to 


«■(  n )  =  0, 


where 

C(n)  =  A -j (n)  -  Xi(n)S^lSo. 

The  matrix  C(n)  is  called  the  Schur  complement  of  9i  in  9 The  equation  (2.3)  has  the  trivial  solution  if  and 
only  if  C(n)  has  full  column  rank.  We  proceed  to  eliminate  the  constraints.  Let 


u(n)  = 


f  u\(n)  \ 

\  wo{n) )  ' 


so  that  ti'i(n)  is  k  x  1  and  Uo(n)  is  p  x  1.  Srnce 

9iU'i(n)  +  S2u-o(n)  =  6, 


we  get 

u'i(n)  =  9f  lb  ~  SilSoWo(n). 
Let 

r(n)  =  -A',(n)9,_l6. 

We  derive 

|[C'(  n ) U'o( n )  -  i  (n)||  =  min. 


(2.4) 


(2.5) 


an  unconstrained  problem  analyzed  in 
complement  matrix  C(n)  recursively  by 


[3],  [4],  Now.  what  about  the  residual  element  rn 


C(n)  = 


(  C(n  -1)\ 

{  r(n)r  /  ' 


Define  the  Schur 
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Partition  the  row  vector  x(n)T  so  that 

r(n)r  =  (  Ji(n)r  12 {n)T  )  . 

where  X[ (n)T  is  1  x  k  and  x 2(n)T  is  1  x  p.  We  get 

c(n)T  =  x2[n)T  -  xi(n)'r5f  lS2. 

Let  v„  denote  the  n-th  element  of  v(n).  The  last  residual  element  of  (2.5)  is  then 

c(n)T  iv2(n)  -  vn  =  x2 (n)Tw2{n)  +  x, {n)Tu'1(n)  +  t>n  -  vn  =  r„, 

i.e.,  the  same  residual  element  as  desired  by  the  constrained  problem  (2.1). 

How  do  we  calculate  r„  recursively?  Suppose  that  we  have  available  a  QR  decomposition  of  the  (n-l)xp 
matrix  C(  n  —  1 ): 

C'(n  -  1)  =  Q(n  -  1  )R(n  -  1). 

where  Q(n  —  1)  is  (n  —  1)  x  p  with  orthonormal  columns  and  the  matrix  R(n  —  1)  is  p  x  p  upper  triangular.  The 
problem  (2.5)  is  reduced  to 


where  u(n  —  1)  =  Q(n  —  1)H  v(n  —  1).  We  triangularize  the  coefficient  matrix  by  a  unitary  matrix  P.  Then 

pH  (  R(n  -  1)  u(n  -  1)\  _  (  R(n)  u(n)\ 

P  V  /“VO7  7  / 

so  that  R(n)  is  p  x  p  upper  triangular.  The  matrix  P  consists  of  p  Givens  matrices.  From  P  and  Q{n  —  1)  we  can 
construct  an  n  x  p  orthonormal  matrix  Q(n)  such  that  C(n)  =  Q(n)R[n)  and  u(n)  =  Q(n)H v(n).  The  desired 
element  rn  is  given  by 

•  •  •  £p  )y . 

where  C\ ,  . .  .  ,cp  denote  cosines  of  the  p  rotations  that  make  up  P. 

3.  PERTURBATION  ANALYSIS 

Elden  [lj  presented  a  perturbation  analysis  of  the  linearly  constrained  least  squares  problem.  Since  his  theory 
is  general,  it  involves  weighted  pseudoinverses  and  their  corresponding  condition  numbers.  In  this  section,  we 
derive  simpler  perturbation  bounds  for  the  solution  w(n)  as  well  as  for  the  residual  element  rn.  To  simplify  our 
presentation,  we  will  drop  the  argument  (n)  for  the  matrices  and  vectors,  and  let  k(M)  denote  the  condition  number 
of  a  matrix  M  with  respect  to  the  2-norm. 

Let  d’  solve  the  perturbed  least  squares  problem 

||  ( ,V|  +  e Ex,  A'z  +  (Ex,  )  u’ll  =  min  (3.1a) 

subject  to  the  perturbed  linear  equality  constraints 

(S\+cEst  S2  +  (E$i)  w  =  b  +  eft,.  (3.1b) 

Suppose  that  t  >  0  is  a  real  variable  and  let 

C  +  tEc  =  (X2  +  tE\2)  —  ( A i  +  tExt  )(S\  -I-  tEs ,  )_l(^2  +  tEs2) 

and 

v  +  tfr  =  -(A',  +tEXlUS]  +tESirl(b  +  tfb). 
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Recall  that  ,St  is  nonsingular  ami  that  C  has  full  column  rank.  Suppose  c  is  sufficiently  small  so  that  for  t  G  [0.( 
we  have  S\  +  t  E-  is  nonsingular  and  C  +  tEr  has  full  column  rank.  Let  u  (t)  solve  the  matrix  equation 


(  S i  ■'S'.1  +  t  j 

v  0  lC  +  tEc)H(C  +  tE') 


w{t)  = 


(  b  +  tfb 

\(C  +  tEc)"(v+  tf,  ) 


(3.2) 


Then  u  ( 0 )  and  tr(f)  are  solutions  to  problems  (2.1)  and  (3.1).  respectively.  Define  w  =  u(0)  and  il  —  u(t).  Then 

(?■  =  ic(0)  -fi(U'(0)  +  0(r). 

Differentiate  (3.2)  with  respect  to  t  and  set  t  =  0.  We  get 

/•:. 


!)’  E»C  + C"  EC  )  ,r(0)  +  (  0  C"t ■ )  <i  (0)  ~  (  EH  v  +  C"f,  ) 


Let 


Then 


.S  [  .Sn 
0  / 


( '  = 


/  0 
0  ( 


.  .  d  = 


and  fd  =  (  ^ 


r"c)  =5-,(r"o-'  and  ||qi<ri|. 


Solving  for  ii(0)  in  (3.3),  we  obtain 


(1(0)  = 


lS  i  ,S  9 

0  c"c 


-1 


^M4r)-(:  4c)-] 


=  S-'(Cr{C)-'C" 


h- 


Es,  Es, 


-  S~l(C"C)-1 


0  Ec 

where  r  =  Cus  —  r  denotes  the  residual  vector.  Furthermore,  by  assuming 

IIMI<  11*11.  II/.  II  <  IMI-  Il^dl  <  IK II 

and 

II  /  r-  r~'  \  II  II  /  L*  C*  \  \\ 

<  \\S\\  IKU 


t-c  r  )  ' 


(Es, 

EsA 

< 

(s'  *:) 

V  0 

Ec) 

V  o  rj 

(3.3) 


(3.4) 


(3.5a) 


(3.5b) 


we  derive  the  inequality 


||«H0)||  <  ||.s-‘||  ||(c//cr‘c"||  [\\d\\  +  IICII  nil!  Hull]  +  p-‘ll  ll(CwC)-1ll  licit  IMI- 

Consequently,  we  obtain  the  following  perturbation  result. 


Lemma.  Using  the  notations  defined  above  and  assuming  that  f  in  (3.1)  is  sufficiently  small  so  that  the  inequalities 
(3.5)  are  satisfied,  we  get 


IK  -  HI 
IMI 


<  < 


|k<S>(C) 


+  k(S')k(02 


[Mj  ] 

IICII  II5II  IM  I!  / 


+  O(C).  Q 


(3.6) 


To  illustrate  the  effect  of  k(S)  on  the  solution  of  (2.1).  consider  a  simple  example  in  which  5  —  (  5(  0  )  and 
.V  =  (  In  In  )■  where  In  is  an  r(  x  n  identity  matrix  and  n  <  h  <  2 n.  By  observation.  «q  =  .Sf‘6  and  u  n  =  -nq. 
Since  k(S)  =  k(S\  )  in  this  example,  we  see  why  the  presence  of  k(.S')  is  necessary  in  (3.6). 
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We  proceed  to  derive  a  bound  for  the  error  in  the  residual.  Let 


m) 


S  i  +  <£5,  So  +  tEs 


j  u'(0  - 


b  +  t  ft, 
f  +  l  ft. 


0  C+tEc 

differentiate  the  equation,  and  then  set  t  =  0.  Using  (3.4)  to  substitute  for  ti'(0),  we  get 


0 

r(0) 


£5,  £s2 

0  £<■• 


=  a-ccf) 


u  + 

V  0 

ESl 

Es, 

0 

Ec 

?)- 


C(CHC) 


Hr\- 1 


£"r 


where  C*  =  (C'HC)  1CH .  Consequently, 

r(0)  =  (/  -  CCt)(Ecw2  -  -  C(CaC)-lEgr. 

As  for  the  residual  element  we  have  rn  =  e^r,  where  en  ~  (0 . 0.1)  denotes  the  n-th  unit  coordinate  vector. 

Using  the  assumptions  (3.5a)  and  noticing  that  us  =  C'^v  and  ;•  =  (/  —  CC^)t’>  we  derive  our  major  result. 


Theorem.  Under  the  same  conditions  as  in  the  Lemma,  we  get 


^  <  <  [||/  -  CC*||(2k(C)  +  D]  +  0(c) 


(3.7) 


and 


f  III7  -  C^IK^C)  +  lien  IIC^II  + 1)]  +  o(c2).  a 


(3.8) 


Here  are  some  additional  remarks.  If  we  set  S  =  7  and  6  =  0,  then  (3.6)  leads  to  a  perturbation  bound  for  the 
standard  least  squares  problem  [2].  We  also  note  that  ||CSu>||2  4-  ||r||2  =  ||d||2.  Thus,  we  can  define 

COS0  =  ||CSuj|/||d|| 

and  use  (1/cosd)  and  tan  9  in  (3.6).  The  bound  (3.7)  is  similar  to  a  result  derived  in  [2].  The  inequality  (3.8) 
indicates  that  |f„  -  rn|  depends  on  «(C)  as  well  as  on  ||r||.  Both  (3.7)  and  (3.8)  can  be  simplified  by  using  the 
relation  that  ||7  —  CC*||  =  min{l,  n  —  p}. 


4.  ERROR  ESTIMATION 

Although  the  error  bound  (3.8)  is  simple,  it  requires  C^en.  whose  computation  involves  at  least  a  back-solve. 
In  this  section,  we  present  an  error  estimation  scheme  for  the  desired  residual  element.  When  the  new  data  vector 
j-(n)r  arrives,  it  is  first  processed  by  S  so  that  Xi(n)T  is  annihilated.  In  particular,  let 


(0) 

1 


4°')  =  x(n)T  and  u'01 


0. 


Then  the  preprocessing  proceeds  as  follows: 


SU 

0 


.(0 

*1+1 


si, 1+1 

ji-il 

•1+ 1 


&l,(J 

Jl-D 


and 


b, 

ii-\) 


\ 

)' 
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for  /  =  1.2 . Jt,  where  <j\  =  "/s/;.  Writing  in  algorithmic  form,  we  have 

for  /  =  1.2 . k 

begin 

i  i-n  / 

<ji  =  ~i  /-v  i  ; 

for  j  =  l  +  1 . q 

Jl)  _  J!-1>  _ 

~i 

=  n</_  1 1  —  gibt 

end. 


The  above  process  shows  that 

(-Jt  +  1  -»*’)  =  -r2(")r  -  J-i(n)T5f152  and  ulk>  =  -x1(n)r5f  lb. 

These  two  variables  are  then  used  for  updating  the  QR  decomposition  of  C(n  —  1)  and  computing  the  residual 
element.  We  present  below  the  algorithm  derived  in  [4], 

for  /  =  1,2 . p 

begin 


(  n )  _  1  Ip  ,  i  lk+l-l)p  , 


a  1  n  -  1 )  ,  ( n  ) 
COS  B,  =  c,  .  /cu  ; 


.  a  |i+(-l)  ,  ( n  I 
sind,  =  zk+l  /c,j  ; 

for  j  =  l  +  1 ,  .  .  p 


begin 


r 


I'D 

lj 


(n-1) 

CIJ 


cos  0/  +  ; 


(*+/- 1) 


sin  9[ : 


jk+n 

~kk] 


(n-l) 

Clj 


sin  0i  +  : 


(k  +  l-l) 
k+j 


cos  0 1 


end; 


i>jn*  =  t'Jn  cos  9\  -h  u*  k+i  Usjn^(. 
u(k  +  l)  =  _[-Jn_1,sin  +  u' *+'->>  cos  01 

end: 


rn  =  «'*+'->n?=lcostfi 


In  the  above,  c(*'  (for  k  —  n  —  l.n)  denotes  the  ( i.  j)-element  of  C(k)  and  rl"'1  the  i-th  element  of  r(k). 


( fc  i 
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Now,  we  discuss  an  error  estimation  scheme  for  the  preprocessing.  Let  '  denote  the  corresponding  computed 
value  and  fl  the  floating  point  computation.  In  the  above  procedure  we  calculate 

9i  =  //(-f“uA./), 
a1'1  =  //(&(,-l)  —  fUgifni))- 

Define  the  relations  between  the  exact  and  computed  quantities  as  follows: 

=  -V  (7 1  ,j  0t  j  (f ) ) . 

gi  =  <//(!+  «i£i(0). 

u‘ '  ’  =  u,,,(l  +  .jl'W',(f)). 
if  =  6/(1  +  ppMO), 

where  |dt,;(c)|  =  0(0,  | ^  *  *  ( € )  |  =  0(c),  |£/(c)|  =  0(c),  |0(,'(O|  =  0(c)  and  |6/(c)|  =  0(c).  The  five  quantities  <r, ;, 
cv/,  r/(i|  and  p/  are  all  real  and  nonnegative.  We  also  assume  that  the  errors  such  as  <Ti,j(j> ij(()  and  t^  fc) 

are  so  small  that  higher  order  terms  like  (c/  /<2>/ /(c))2  and  (<x//0//(c ))(<,'{  —  1 ' ir j ^  "(«))  can  he  ignored.  Using  the 
lemma  in  [3],  we  obtain  the  following  algorithm  for  estimating  the  errors  in  preprocessing. 

for  /  =  1,2,  ....  k 


«/=max{Cjl  1,.(T,/}; 

f-  ;  -  •' 1 1 . 4 


,|()  I--'1  1  max(«i.<»i., ) 


//)  _  |u“  "'/''  ''|  +  |^i»[  max)q,  /i,} 

'<  uTm 


As  explained  in  [3],  the  above  estimation  scheme  can  be  incorporated  with  the  preprocessing  procedure  and  imple¬ 
mented  on  the  same  systolic  architecture.  Additional  time  is  minimal  because  the  calculations  can  be  carried  out 
during  the  otherwise  idle  time  of  the  processors. 

The  error  estimate  for  rn  can  be  obtained  by  the  algorithm  presented  in  [3j  using  (C/t+'i  '  '  <."«*')  and  q{k)  as 

the  error  estimates  for  (i^A  and  u{k),  respectively.  Again,  we  list  the  error  estimation  algorithm  and 

refer  the  details  to  [3],  Define  the  relations  between  the  exact  and  computed  quantities  as  follows: 


i)*’  =  ^‘’(l+C/V^’tc)). 

fi“>  =  u'^l  +V*’fl“,(c)). 

- ( n  —  i  j  =  +^l  p+1a1.f,+  l(c)). 


c["’  =  c‘"’(l  +  at]0t](c)). 

COS0|  =  COS0|(1  +  CT/./^A/Ac))* 

sin  §i  =  sin  0/(1  +  <r( ./^’"/(c)), 

,<n|  (n|.  , 

(1  +  <T,  r+{0,  ,p+i(c)). 


rn  =  rn(l  +  T7<?(e ) ) 
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The  following  algorithm  estimates  the  error  in  the  last  element  of  the  residual  vector: 
for  /  =  1.2 . v 


I  r  .{k  +  l-  1  j  1 

/  =  n\ax{fiu.isk+l  j 
for  j  =  l  +  1 .  .  .  .  p 


I  *’  »»  —  l  )  A  Is-  1  I  |  I  k  +  f  ~ *  l  t  rt  |  .  |  +1  —  1  )  I  I 

_  ly,  ,  -n  sin  max  { ^  +  /  '’nil 

1  1  'os  +- j*++|,_  "  sln#'1 

.  i  n  -  n  j  j  /  till  (  k  + 1  —  I  i  .  I  •  1  k  +1  —  |  j 

.(/+  k  )  _  \ri  J  Sin*,  inaxK!  }I  +  |;k+i  m:»«  Uk+  y  Ti  . 


_  1 1  '  ”  '  “  cos  H,  Iiu»|(i  p  +  l  Ot, I  H  +  lu1*  *'  "  sin  1*1  max(  V*4'  -  11 

~ *  |t  j  1  ’  cos  (*i+ul‘+,~1'sin»i| 

(  k  +  ( )  _  I* ”  «in  *i  rtiaxlfi  ,  +  i  -n.i  } I -*- 1 u* 1,41  ~  11  ‘-os  »,  manj  r ?  1  <■  -M  -  1  >  •r,  ( jj 

^  It4*”1  sin  <*i -u|l,+*~  u  <~os  <*t| 


0  =  r;“'*-rl  niax{(T,  1  ,  ,(Tpp}  . 

5.  AN  EXAMPLE 

The  example  in  this  section  shows  that  the  computed  residual  element  may  he  accurate  even  when  the  matrix 
C  is  ill-conditioned.  In  this  case,  the  proposed  scheme  gives  a  better  error  estimate  than  (3.8).  Both  the  MS 
algorithm  and  tne  error  estimation  algorithm  were  implemented  using  MATLAB  and  run  on  a  VAX  8300  with 
machine  precision  f  =  1  1102  a  10”lh  in  the  Communications  Research  Laboratory  at  McMaster  University. 


Example.  Suppose  the  exact  constraint  matrix  and  corresponding  right  side  vector  are 

(  103  0  0  1  0  0\  / -1200UUv2/"\ 

0  10-3  0  0  1  0  and  b  =  v'10/700  J 

0  0  1  0  0  1  /  \  6v/5/7  / 


Thus  we  set  the  error  estimates  as  <r,  j  =  /(/  =  1  and  <*>,  j ( e )  =  t'i(f)  =  (.  The  data  matrix  at  time  n  —  1  is 

/  — 1  -n/5  — 2  vTo  0  0  0\ 

A  ( n  -  1)  =  [  0  -1  y/2  000 


0  0 


-1  0  0  0 


Suppose  we  know  the  exact  R(n  —  1)  and  u(n  —  1): 

(  0.001  1000\/5  2vTo  \  /  — 10 >/2 /7 

0  1000  -2v/2  and  u(n  —  1)  =  J  4vTo/7 

001/  \  61/5/7 

Similarly,  the  error  estimates  of  their  elements  are  all  initialized  as  t.  Now  the  new  data 


j[(n)T  =  (-l  -y/E  -2v/10  0.001  0  0)  and  «((”  =  0 
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ar<  available  and  their  error  estimates  are  initialized  as  ^  =1.  for  7  =  1.  .  6.  and  if  =  1.  respect iv  -ly.  After 

preprocessing.  we  yet  r(n)r  =  (  0  002  lOOOv'5  2v/10  )  and  >  „  =  -10  The  corresponding  error  estimates 

are  =  1.  for  j  =  4.  5.6.  and  if''  =  23.  The  QR  updating  scheme  and  its  error  estimation  algorithm  are  then 
applied  to  Rln  —  1).  u(n  -1).  r(»i).  i and  their  error  estimates.  The  exact  residual  element  r,:  =  Cv  2  35.  The 
computed  error  is 

jr„  -  r„|  =  1.11  s  10-i,; 

The  condition  number  of  (  '(  u )  is  4  6  <  10'’ and  the  error  bound  as  given  by  (3.8)  equals  3  40  ■  10“'.  The  estimation 
algorithm  gives  a  much  more  accurate  value  of  9  62  «  10“  |f  . 
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